"""
    test_statevectors.jl

Test suite for the StateVectors module of GSICoreAnalysis.jl.
This module tests state vector operations including allocation,
variable access, transformations, and I/O operations.
"""

using Test
using GSICoreAnalysis
using GSICoreAnalysis.StateVectors
using LinearAlgebra
using Random

# Set reproducible seed
Random.seed!(1011)

@testset "StateVectors Module Tests" begin
    
    @testset "StateVector Construction" begin
        # Test basic construction
        config = AnalysisConfig{Float64}(grid_size = (6, 4, 3))
        sv = StateVector(config)
        
        @test sv isa StateVector{Float64}
        @test sv.config === config
        @test sv.is_allocated == true
        @test length(sv.values) > 0
        @test sv.grid.nx == 6
        @test sv.grid.ny == 4
        @test sv.grid.nz == 3
        @test sv.grid.total == 6 * 4 * 3
        
        # Test precision specification
        config32 = AnalysisConfig{Float32}(grid_size = (4, 4, 2))
        sv32 = StateVector{Float32}(config32)
        @test sv32 isa StateVector{Float32}
        @test eltype(sv32.values) == Float32
        
        # Test different grid sizes
        small_config = AnalysisConfig{Float64}(grid_size = (2, 2, 1))
        large_config = AnalysisConfig{Float64}(grid_size = (12, 8, 4))
        
        sv_small = StateVector(small_config)
        sv_large = StateVector(large_config)
        
        @test length(sv_small.values) < length(sv_large.values)
        @test sv_small.is_allocated == true
        @test sv_large.is_allocated == true
    end
    
    @testset "Memory Management" begin
        config = AnalysisConfig{Float64}(grid_size = (4, 4, 2))
        sv = StateVector(config)
        
        # Test initial allocation
        @test sv.is_allocated == true
        original_length = length(sv.values)
        @test original_length > 0
        
        # Test variable views are properly sized
        @test size(sv.u) == (4, 4, 2)
        @test size(sv.v) == (4, 4, 2)
        @test size(sv.t) == (4, 4, 2)
        @test size(sv.q) == (4, 4, 2)
        @test size(sv.oz) == (4, 4, 2)
        @test size(sv.cw) == (4, 4, 2)
        @test size(sv.ps) == (4, 4)
        @test size(sv.stl) == (4, 4)
        @test size(sv.sti) == (4, 4)
        
        # Test deallocation
        deallocate_state(sv)
        @test sv.is_allocated == false
        @test length(sv.values) == 0
        
        # Test re-allocation
        allocate_state(sv)
        @test sv.is_allocated == true
        @test length(sv.values) == original_length
        
        # Test double allocation (should be safe)
        allocate_state(sv)
        @test sv.is_allocated == true
        @test length(sv.values) == original_length
        
        # Test double deallocation (should be safe)
        deallocate_state(sv)
        deallocate_state(sv)
        @test sv.is_allocated == false
    end
    
    @testset "Variable Access and Views" begin
        config = AnalysisConfig{Float64}(grid_size = (3, 3, 2))
        sv = StateVector(config)
        
        # Test that views are properly connected to main values array
        original_values = copy(sv.values)
        
        # Modify through variable view
        sv.u[1, 1, 1] = 99.5
        @test sv.values != original_values  # Main array should be modified
        
        # Modify through main values array
        sv.values[1] = -77.3
        @test sv.u[1, 1, 1] == -77.3  # View should reflect the change
        
        # Test all variable views\n        var_names = [\"u\", \"v\", \"t\", \"q\", \"oz\", \"cw\", \"ps\", \"stl\", \"sti\"]\n        \n        for var_name in var_names\n            var_slice = get_variable_slice(sv, var_name)\n            @test var_slice isa AbstractArray\n            @test eltype(var_slice) == Float64\n            \n            if var_name in [\"u\", \"v\", \"t\", \"q\", \"oz\", \"cw\"]\n                @test ndims(var_slice) == 3\n                @test size(var_slice) == (3, 3, 2)\n            else  # 2D variables\n                @test ndims(var_slice) == 2\n                @test size(var_slice) == (3, 3)\n            end\n        end\n        \n        # Test invalid variable name\n        @test_throws Exception get_variable_slice(sv, \"invalid_variable\")\n    end\n    \n    @testset \"Variable Setting and Getting\" begin\n        config = AnalysisConfig{Float64}(grid_size = (3, 2, 2))\n        sv = StateVector(config)\n        \n        # Test setting 3D variable\n        test_temp = ones(3, 2, 2) * 285.0\n        set_variable_slice!(sv, \"t\", test_temp)\n        retrieved_temp = get_variable_slice(sv, \"t\")\n        @test retrieved_temp ≈ test_temp atol=1e-12\n        \n        # Test setting 2D variable\n        test_ps = ones(3, 2) * 101325.0\n        set_variable_slice!(sv, \"ps\", test_ps)\n        retrieved_ps = get_variable_slice(sv, \"ps\")\n        @test retrieved_ps ≈ test_ps atol=1e-12\n        \n        # Test dimension mismatch\n        wrong_size_data = ones(2, 2, 2)  # Wrong size for this grid\n        @test_throws AssertionError set_variable_slice!(sv, \"t\", wrong_size_data)\n        \n        # Test get_state_pointer (alias for get_variable_slice)\n        temp_ptr = get_state_pointer(sv, \"t\")\n        @test temp_ptr === get_variable_slice(sv, \"t\")\n    end\n    \n    @testset \"State Vector Operations\" begin\n        config = AnalysisConfig{Float64}(grid_size = (4, 3, 2))\n        sv1 = StateVector(config)\n        sv2 = StateVector(config)\n        sv3 = StateVector(config)\n        \n        # Fill with test data\n        randn!(sv1.values)\n        randn!(sv2.values)\n        \n        # Test dot product\n        dot1 = dot_product_state(sv1, sv2)\n        dot2 = dot_product_state(sv2, sv1)  # Should be symmetric\n        \n        @test dot1 isa Float64\n        @test dot1 ≈ dot2 atol=1e-12  # Symmetry\n        \n        # Test with self\n        self_dot = dot_product_state(sv1, sv1)\n        @test self_dot >= 0.0  # Positive semi-definite\n        @test self_dot ≈ norm_state(sv1)^2 atol=1e-10\n        \n        # Test AXPY-style operation: sv2 = sv2 + α*sv1\n        α = 2.5\n        sv2_orig = StateVector(config)\n        assign_state!(sv2_orig, sv2)\n        \n        add_state!(α, sv1, sv2)\n        \n        # Verify the operation\n        expected = sv2_orig.values + α * sv1.values\n        @test sv2.values ≈ expected atol=1e-12\n        \n        # Test assignment operation\n        assign_state!(sv3, sv1)\n        @test sv3.values ≈ sv1.values atol=1e-15\n        \n        # Modify sv1 and ensure sv3 is independent\n        original_sv3_vals = copy(sv3.values)\n        sv1.values .+= 1.0\n        @test sv3.values ≈ original_sv3_vals atol=1e-15\n        \n        # Test copy_state! (alias for assign_state!)\n        copy_state!(sv3, sv2)\n        @test sv3.values ≈ sv2.values atol=1e-15\n    end\n    \n    @testset \"Zero and Norm Operations\" begin\n        config = AnalysisConfig{Float64}(grid_size = (3, 3, 2))\n        sv = StateVector(config)\n        \n        # Fill with non-zero values\n        randn!(sv.values)\n        @test !all(sv.values .== 0.0)\n        \n        # Test zero operation\n        zero_state!(sv)\n        @test all(sv.values .== 0.0)\n        \n        # Test norm with known values\n        fill!(sv.values, 2.0)\n        norm_val = norm_state(sv)\n        expected_norm = sqrt(length(sv.values) * 4.0)  # 2.0^2 * n_elements\n        @test norm_val ≈ expected_norm atol=1e-10\n        \n        # Test norm of zero vector\n        zero_state!(sv)\n        @test norm_state(sv) ≈ 0.0 atol=1e-15\n    end\n    \n    @testset \"Variable-Specific Operations\" begin\n        config = AnalysisConfig{Float64}(grid_size = (4, 3, 2))\n        sv = StateVector(config)\n        \n        # Set up realistic atmospheric values\n        nx, ny, nz = 4, 3, 2\n        \n        # Temperature field (decreasing with height)\n        for k in 1:nz, j in 1:ny, i in 1:nx\n            height_factor = (k - 1) / (nz - 1)\n            sv.t[i, j, k] = 288.0 - 40.0 * height_factor\n        end\n        \n        @test all(sv.t[:, :, 1] .> sv.t[:, :, 2])  # Surface warmer than upper level\n        @test minimum(sv.t) >= 240.0  # Reasonable temperature range\n        @test maximum(sv.t) <= 300.0\n        \n        # Wind field\n        sv.u .= randn(nx, ny, nz) * 10.0  # ±10 m/s winds\n        sv.v .= randn(nx, ny, nz) * 5.0   # Smaller v component\n        \n        @test abs(maximum(sv.u)) <= 50.0  # Reasonable wind speeds\n        @test abs(maximum(sv.v)) <= 25.0\n        \n        # Humidity (positive values only)\n        sv.q .= abs.(randn(nx, ny, nz)) * 0.01\n        @test all(sv.q .>= 0.0)  # Humidity should be non-negative\n        @test maximum(sv.q) <= 0.05  # Reasonable humidity values\n        \n        # Surface pressure\n        sv.ps .= 101325.0 .+ randn(nx, ny) * 1000.0\n        @test all(sv.ps .> 80000.0)   # Reasonable pressure range\n        @test all(sv.ps .< 110000.0)\n        \n        # Other variables\n        fill!(sv.oz, 3e-6)     # Typical ozone mixing ratio\n        fill!(sv.cw, 0.0)      # No cloud water\n        fill!(sv.stl, 285.0)   # Land surface temperature\n        fill!(sv.sti, 271.0)   # Ice surface temperature\n        \n        @test all(sv.oz .> 0.0)\n        @test all(sv.cw .== 0.0)\n        @test all(sv.stl .> 250.0)\n        @test all(sv.sti .> 250.0)\n    end\n    \n    @testset \"File I/O Operations\" begin\n        config = AnalysisConfig{Float64}(grid_size = (3, 2, 2))\n        sv_orig = StateVector(config)\n        \n        # Fill with distinctive values for each variable\n        nx, ny, nz = 3, 2, 2\n        \n        # Create distinctive patterns for each variable\n        for k in 1:nz, j in 1:ny, i in 1:nx\n            sv_orig.u[i, j, k] = 10.0 + i + j + k\n            sv_orig.v[i, j, k] = 20.0 + i - j + k\n            sv_orig.t[i, j, k] = 280.0 + i * j * k\n            sv_orig.q[i, j, k] = 0.01 * (i + j + k)\n            sv_orig.oz[i, j, k] = 3e-6 * i * j\n            sv_orig.cw[i, j, k] = 0.001 * k\n        end\n        \n        for j in 1:ny, i in 1:nx\n            sv_orig.ps[i, j] = 101000.0 + i * j * 100.0\n            sv_orig.stl[i, j] = 285.0 + i - j\n            sv_orig.sti[i, j] = 271.0 + j - i\n        end\n        \n        # Write to file\n        test_file = tempname()\n        write_state(sv_orig, test_file)\n        \n        @test isfile(test_file)\n        \n        # Check file content has reasonable structure\n        file_content = read(test_file, String)\n        @test contains(file_content, \"# GSI State Vector\")\n        @test contains(file_content, \"# Variable: u (3D)\")\n        @test contains(file_content, \"# Variable: ps (2D)\")\n        \n        # Read back\n        sv_read = StateVector(config)\n        read_state(sv_read, test_file)\n        \n        @test sv_read.is_allocated == true\n        \n        # Note: The read_state implementation is simplified and may not\n        # perfectly reconstruct all values, but it should work for basic cases\n        \n        # Clean up\n        rm(test_file)\n    end\n    \n    @testset \"Grid Dimensions\" begin\n        # Test GridDimensions construction\n        grid1 = GridDimensions(10, 8, 5)\n        @test grid1.nx == 10\n        @test grid1.ny == 8\n        @test grid1.nz == 5\n        @test grid1.total == 10 * 8 * 5\n        \n        # Test construction from tuple\n        grid2 = GridDimensions((6, 4, 3))\n        @test grid2.nx == 6\n        @test grid2.ny == 4\n        @test grid2.nz == 3\n        @test grid2.total == 6 * 4 * 3\n        \n        # Test in StateVector context\n        config = AnalysisConfig{Float64}(grid_size = (5, 4, 3))\n        sv = StateVector(config)\n        \n        @test sv.grid.nx == 5\n        @test sv.grid.ny == 4\n        @test sv.grid.nz == 3\n        @test sv.grid.total == 5 * 4 * 3\n    end\n    \n    @testset \"Error Handling\" begin\n        config = AnalysisConfig{Float64}(grid_size = (3, 3, 2))\n        sv1 = StateVector(config)\n        sv2 = StateVector(config)\n        \n        # Test operations on unallocated vectors\n        deallocate_state(sv1)\n        deallocate_state(sv2)\n        \n        @test_throws AssertionError dot_product_state(sv1, sv2)\n        @test_throws AssertionError add_state!(1.0, sv1, sv2)\n        @test_throws AssertionError assign_state!(sv1, sv2)\n        @test_throws AssertionError norm_state(sv1)\n        @test_throws AssertionError zero_state!(sv1)\n        \n        # Test I/O on unallocated vector\n        @test_throws AssertionError write_state(sv1, \"dummy.txt\")\n        @test_throws AssertionError read_state(sv1, \"nonexistent.txt\")\n        \n        # Test variable access on unallocated vector\n        @test_throws AssertionError get_variable_slice(sv1, \"t\")\n        @test_throws AssertionError set_variable_slice!(sv1, \"t\", ones(3, 3, 2))\n    end\n    \n    @testset \"Numerical Stability\" begin\n        config = AnalysisConfig{Float64}(grid_size = (4, 4, 2))\n        \n        # Test with very small numbers\n        sv_small = StateVector(config)\n        fill!(sv_small.values, 1e-15)\n        \n        norm_small = norm_state(sv_small)\n        @test norm_small >= 0.0\n        @test isfinite(norm_small)\n        \n        # Test with very large numbers\n        sv_large = StateVector(config)\n        fill!(sv_large.values, 1e10)\n        \n        norm_large = norm_state(sv_large)\n        @test isfinite(norm_large)\n        \n        # Test dot product doesn't overflow\n        dot_val = dot_product_state(sv_large, sv_large)\n        @test isfinite(dot_val)\n        @test dot_val >= 0.0\n        \n        # Test mixed operations\n        add_state!(1e-10, sv_small, sv_large)\n        @test all(isfinite.(sv_large.values))\n    end\n    \n    @testset \"Performance Considerations\" begin\n        # Test that operations are reasonably efficient\n        config = AnalysisConfig{Float64}(grid_size = (20, 15, 8))\n        sv1 = StateVector(config)\n        sv2 = StateVector(config)\n        \n        randn!(sv1.values)\n        randn!(sv2.values)\n        \n        # These operations should complete quickly\n        @test (@elapsed dot_product_state(sv1, sv2)) < 1.0  # Less than 1 second\n        @test (@elapsed add_state!(1.5, sv1, sv2)) < 1.0\n        @test (@elapsed assign_state!(sv2, sv1)) < 1.0\n        @test (@elapsed norm_state(sv1)) < 1.0\n        @test (@elapsed zero_state!(sv2)) < 1.0\n        \n        # Variable access should be efficient\n        @test (@elapsed get_variable_slice(sv1, \"t\")) < 0.1\n        @test (@elapsed set_variable_slice!(sv1, \"u\", sv1.u)) < 0.1\n    end\n    \n    @testset \"Memory Layout and Views\" begin\n        config = AnalysisConfig{Float64}(grid_size = (3, 2, 2))\n        sv = StateVector(config)\n        \n        # Test that variable views are actually views of the main array\n        original_ptr = pointer(sv.values)\n        \n        # The views should point to parts of the main array\n        u_ptr = pointer(sv.u)\n        v_ptr = pointer(sv.v)\n        \n        @test u_ptr >= original_ptr\n        @test v_ptr > u_ptr  # v should come after u in memory\n        \n        # Test that modifying views affects main array\n        original_checksum = sum(sv.values)\n        sv.t[1, 1, 1] += 100.0\n        new_checksum = sum(sv.values)\n        \n        @test abs(new_checksum - original_checksum - 100.0) < 1e-12\n        \n        # Test variable ordering in memory\n        # Expected order: u, v, t, q, oz, cw, ps, stl, sti\n        fill!(sv.values, 0.0)\n        \n        sv.u[1, 1, 1] = 1.0\n        sv.v[1, 1, 1] = 2.0\n        sv.t[1, 1, 1] = 3.0\n        \n        # Find where these values appear in the main array\n        u_idx = findfirst(x -> x == 1.0, sv.values)\n        v_idx = findfirst(x -> x == 2.0, sv.values)\n        t_idx = findfirst(x -> x == 3.0, sv.values)\n        \n        @test u_idx < v_idx < t_idx  # Verify ordering\n    end\nend\n\nprintln(\"StateVectors module tests completed successfully.\")"